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SUMMARY 

An analytic solution for the critical, monoenergetic , bare, infinite 
cylinder is presented. The solution is obtained by modifying a previous 
development based on a neutron density transform and Case's singular 
eigenfunction method. Numerical results for critical radii and the neu- 
tron density as a function of position are included and compared with the 
results of other methods. 


I. INTRODUCTION 

The exact solutions for the critical dimensions of multiplying media 
and the associated neutron density distributions can be used as standards 
in evaluating the more widely used approximate methods of transport 
analysis such as discrete ordinates or spherical harmonics. Of interest 
here is the exact solution for the monoenergetic, bare, infinite cylinder. 
Critical radii for monoenergetic bare cylinders have been determined by 
several methods. Carlson and Bell Q.) interpolated between values calcu- 
lated with the extrapolated endpoint and variational methods to derive 
what have long served as reference values. More recently, approximate 
solutions have been obtained by Hendry, (2) using Fourier expansions, and 
by Hembd,Q) using a Fourier transformation of the integral equation 
(I Tn Method). Mitsis, using neutron density transforms, was able to 
apply Case's method (-D to obtain exact analytical solutions for spheres 
and cylinders. Gibbs generalized this work to an arbitrary convex 
geometry and showed a reduction to the Mitsis solutions for the two 
specialized geometries. Numerical results for critical spheres were pre- 
sented by Mitsis; however, a recent bibliography (Z) of transport theory 
does not include any reference to numerical results derived from the ana- 
lytical solution for the critical cylinder. Initial investigation under 
the present study indicated that the Mitsis formulation for critical 
cylinders is not convergent. 

The solution of the critical cylinders presented here was obtained 
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by reformulating the spatial function and the continuum expansion coef- 
ficient of the Mitsis solution to achieve convergence. Numerical re- 
sults for critical radii and the asymptotic and total neutron densities 
as a function of position are included. 


II. ANALYSIS 

This analysis includes the development of a neutron density trans- 
form for the infinite, rotationally invariant cylinder, the solution of 
the transformed transport equation by the Case singular eigenfunction 
expansion technique, and the application of this solution to the criti- 
cal cylinder problem. The analysis follows that of Mitsis, the primary 
features of which are included here for completeness and to provide a 
background for introducing the modifications necessary for obtaining a 
solution. 


A. Development of Neutron Density Transform 


In general, the spatial distribution of the neutron density is 
given by Eq. (1) where c denotes the mean number of neutron second- 
aries per collision. 


P (r) 



p(^)e-' £ - £ ' A , 


r - r 


i 1 2 


( 1 ) 


Considering the geometry of interest shown in Fig. 1, if we place the 
position vector x_ along the x-axis and represent r_' by the cylin- 
drical coordinates (t, ct, a), we can obtain the following relationships: 


2 2 2 

x = r + t - 2rt cos a 


x + z 


= !r - r'|2 


|r - r’ ) 2 = r 2 + t 2 - 2rt cos a + z 2 


(2a) 

(2b) 

(2c) 


Substitution of Eq. (2c) into Eq. (1) yields 


P(r) 


pR 2 tt 


2 "J, 


tp (t)dt 



da 


exp 


[-/r 2 + t 2 - 2rt 


cos a + z 


dz 


2 2 2 
r 4- t - 2rt cos a + z 


( 3 ) 


The integral over z is reduced in the following manner: 
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r\ 

exp! 

/ 2 2 0 „ , 2 
lr/r + t - 2rt cos a + z 

dz 

r 

Jl . 2 . 

^-vx + z Jdz 

1 

2 2 2 1 
r + t 2rt cos a + z ^ 

2,2 

0 x + 2 


(4a) 


>CrxA + ( 


z/x) 2 J 


1 / exp 

x 2 vJ n 1 + (z/x) 2 


dz 


(4b) 



(4c) 


Changing the order of integration, letting v = z/x, and using an iden- 
tity^) for K Q (u) , the modified Bessel function of the second kind, 

Eq. (4c) becomes 


1_ 

x 


O' 

du / SH 

'0 


C-u/l + V 2 


' x v — fi + v 2 

Finally, letting y = u/x, Eq. (4d) becomes 


dv = - 
x 



K (u)du 
o 


(4d) 


K o(^ 


2 2 
r + t - 2rt cos a jdy 


Substitution of Eq. (4e) into Eq. (3) gives 

\R /^,2it 


(4e) 


>(r) = tp(t)dt da + t 2 - 2rt cos a^dy 


(5) 


The integral over a is performed by first splitting the integral over 
t into two integrals ranging from 0 to r and from r to R and 
then applying the addition theorem(_) for K 0 . 



2rt cos 


OO 

• -2 


K (r) I (t) 
ina ) n n 

K (t)I (r) 
n n 


r >_ t 
r < t 


(6) 


The exponential in Eq. (6) leads to sine and cosine terms for all orders 
other than zero. These terms integrate to zero in Eq. (5) and we are 
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left with a 2 tt contribution for the integral over a . Thus p (r) 
becomes 


P (r) = c 


s' *rr 


K Q (ry) 1 Q ( ty ) tp(t)dt 



K o (ty)I o (ry)tp(t)dt 


( 7 ) 


Making the substitution y = .l/y we obtain 


P (r) 



K Q (r/y) I Q (t/y) tp(t)dt 


K o (t/y)I o (r/y)tp(t)dt 


( 8 ) 


Eq. (8) is Eq. (5.5-1) of Mitsis.^A) It is in the form of a kernel 
K(r,t) which is related to a Green's function G(r,t;y) 


P (r) = c 



tp (t)K(r ,t)dt 


0 < t < R 


(9) 


K(r,t) 



~ G(r,t;y)dy 
y 


G(r,t;y) = 


K (r/y)I (t/y) 
o o 


K (t/y)I (r/y) 
O o 


r > t 
r < t 


(9 a) 


(9b) 


which is associated (.12.) as indicated below with the operator 
[d^/dr^ + l/r(d/dr) - l/y^]. Next Mitsis introduces a pseudo neutron 
distribution function $(r,y) by 


$(r,y) 



tp (t)G(r ,t ;y)dt 


( 10 ) 



tp(t)K o (r/y)I Q (t/y)dt 


+ c 


tp(t)K o (t/y)I Q (r/y)dt 


(10a) 


o 

Operation on $(r,y) by multiplying (1/y ) and integration over y re- 
turns the neutron density. 



5 


Differentiation 
relation (§) 



of Eq. (10a) with respect to r 

K n + l (z) V 2) + W z) V z) 


and use of the 
= 1/z 


( 11 ) 


yields 


3$(r ,y) 1 3$(r ,u) 

a 2 r 3r 

3r 


V 


r 1 $(r ,y ' )dy 1 

Jo u’ 2 


( 12 ) 


The boundary conditions for Eq. (12) follow from the definition for 
3>(r,y) in Eq. (10a). Setting r = 0 causes the first integral to 
vanish as the second integral remains finite. At the outer boundary 
we have 


K 0 (R/y)*'(R,y) + ^ K 1 (R/y)*(R,y) - 0 y>0 (13) 

Eqs. (10), (11) define a transform pair with the function 3>(r,y) obeying 
the inhomogeneous differential equation (12) . The inhomogeneity in 
Eq. (12) is proportional to the neutron density by Eq. (11) . 


B. Development of Eigenfunctions by Case's Method 

Mitsis goes on to determine a complete set of eigenfunctions for 
this transform. He assumes a set of separable solutions of the form: 


$ v (r,y) = R v (r)n v (y) (14) 

Substitution of Eq. (14) into Eq. (12) yields 

R^'(r) + i R^(r) - ~ R (r) = 0 (15) 

v 

and 

- ^ 2 ^n v (y) = hyCp^dy'; y > 0 (16) 

o 

with (l/v^) being the separation constant. With the boundary condition 
of R v (o) finite, the solution of Eq. (15) is restricted to the modified 
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Bessel function of the first kind. 

R v (r) = I 0 (r/v) (17) 

Following Case's method, n v (y) is normalized by taking 

\ n v (y)dy = 1 (18) 

y 

after which Eq. (16) is written as 

(v 2 - y 2 )n y ()j) = cv 2 y 2 ; 0 < y < 1 (19) 

For v in the interval (-1,1), singularities occur at v = ±y . Account- 
ing for these singularities, the solution of Eq. (19) becomes 

p 2 2 

n v (»)::=' 2 2 + 5<v)«(v - V) + 5(V)6(V + y) (20) 

v - y 

where P indicates that principal values are to be taken when integrat- 
ing expressions containing n v (y) and £(v),£(v) are arbitrary functions. 
It can be seen in Eq. (19) that P v (y) = n v (-y). Applying this symmetry 
in Eq. (20) requires that ?(v) = £(v) . Thus we can. replace these func- 
tions in Eq. (20) by v^-X(v) to obtain 

2 2 

n v (y) = C , 2 V y 2 + y2 ^(v)jj(y - v) + <5(y + v)] ; y > o (21) 
v - y 

For v not in (-1,1) the solution of Eq. (19) is 


( 22 ) 


(22a) 

The discrete and continuum eigenfunctions ng(y) and n v (y) are directly 
related to those developed by Case for plane geometry, ^ + (y) and 
<t> v (y), and share their completeness properties. 

n Q (y) = y 2 t4> 0+ (y) + 4> 0 _(y) ] 


2 2 

f?V 

V y) = ~T — s 

v 0 " y 


where ±Vq are the roots of 


cv. tanh "'"(l/v) = 1 



(23a) 
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n v (p) = y 2 [4> v (p) + <L v (y)] 


where 


and 


, / s _ £ 0 

*0± M - 2 v Q |y 


P 

<t> v (y) = § ~ Z ~ + X(v)6(y - v) 


(23b) 


(23c) 


(23d) 


A(v) = 1 - cv tanh 


(23e) 


Thus an expansion for an even function $(y) defined in the range 
-1 < y < 1: 


$(y) = a Q n 0 (y) + 



A(v)n y (y)dv 


= Vo (y) + 



A(v)n v (y)dv; 


n -v (y):= n v (y) 

A(v) = A(v) + A(-v) 


(24) 


can be rewritten as: 



$(y) = y <,a Q [(j> 0+ (y) + <j> 0 _(y)] + J A(v) [cj> v (y) + cj>_^ (y ) ] dv > (25) 


Eq. (25) is an expansion of the even function (y) = 3>(y)/y in terms 
of the complete set <j> v (y) with a unique set of coefficients ag and 
A(v). We now insert Eqs. (24) and (17) into Eq. (14) to form the gen- 
eral solution of Eq. (12). 


$(r,y) = a 0 n 0 (y)I 0 (r/v 0 ) + J A(v)n v (y)Ig(r/v)dv 


(26) 


Operation on $(r,y) as in Eqs. (11) and (18) returns the neutron 
density. 


p(r) 


a oV r/ V 



A(v)l Q (r/v)dv 


(27) 
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C. Application to the Critical Cylinder 

To obtain the solution for the critical cylinder we evaluate 
<J>(r,y) and its first derivative at r = R and insert them into the 
boundary condition, Eq. (13), obtaining: 



a 0 n 0 (u)q(v 0 ,y) + j A(v)n v (y)q(v,y)dv =0; y _> 0 (28) 


where 


q(v,y) = R 


i K Q (R/y)I 1 (R/v) + i K^R/yOl^R/v) 


(29) 


The general procedure is to separate Eq. (28) into a singular part which 
can be reduced to a Fredholm type of integral equation to determine the 
expansion coefficient and a corresponding nonsingular expression for the 
eignefunction. Toward this end Mitsis introduces the function 


H(v,y) 


_ q(v,y) - 1 


v - y 


(30) 


which was found by the authors to give a diverging solution. Here we 
modify his procedure by reformulating H(v,y) as 


H'(v,y) = X 


"I 0 (R/y) 


v - y (I 0 (R/v) 


q(v,y) - 1 


(31) 


1 C R MR/v) R ^ 

- y^>y*Ao t^F) + y io (R/ ^ )K i (R/ ^ - j (31a) 

Considering the behavior of the modified Bessel functions it can be seen 
that H(v,y) is an unbounded function for increasing R/ v whereas 
H (v,y) is bounded for all variables. 


y R/ v) r A 

q (v,y) = ~ f ' (R/y) V 1 + ( v - y)H'(v,y)J 


(32) 


We also redefine the expansion coefficient. 


A (v) = A(v) I q (R/v) 


(33) 


Substitution of Eqs. (32) and (33) into Eq. (28) along with the use of 
Eqs. (23 a-d) completes the desired separation. 
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A* ( v) 4> ^ (y )<i v = $ (y); y >_ 0 (34) 

$'(y) is the nonsingular part defined by 

3>'(y) = -a Q [<j) 0 + (y) + <j) 0 _(y)]l 0 (R/y)q(v 0 ,y) 

1 £ vA* ( v)C2vh' ( v,y) + l] , 

2 V + y 

(35) 

Reduction of (34) to a Fredholm equation proceeds as in the solution for 
critical s lab s » — ?.) and results in 




a' ( y) = X(y)g(c,y)$ ' (y) 


x~(y)X (y) 



y(v)$ (v)dv 
v - y 


(36) 


and the auxiliary condition 



y(y)$ (y)dy = 0 


(37) 


which is the criticality condition for this geometry. Additional func- 
tions and their boundary values introduced by this reduction are standard 
forms (A 2 ) with the following definitions. 

g(c,y) = 1 /{x 2 (y) + (c 2 TT 2 y 2 /4)} (38a) 

1/x (y)X + (y) = g(c,y)(v 2 - y 2 ) (1 - c)x(-y) (38b) 

y(y)=f ii — 5 (38c) 

(v Q - y ) (1 - c)x(-y) 

x(-z) = K(-z) /{l/x(0) 4- z} (39a) 


K(-z) = 1 - 



[1 - x 2 ( 0 )y 2 ]dy 


1 — 2 )K( - y) (y + z) 


(39b) 
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x 2 (0) = 1/ v g (1 ' c) 


(39c) 


Eqs. (35), (36), and (37) form the solution. Successive approximations 
are constructed by the following iteration scheme: 


Aq(p) = 0 


$q(p) = -a 0 (j> 0+ (y) + <j) 0 _(y)]l 0 (R/y)q(v 0 ,y) 

, Y(v)4>'(v)dv 

A (y) = l(y)g(c,y)$ (y) / 

n+1 n x-(y)X + (y)vi) 


(40a) 

(40b) 

(40c) 


n 



(40d) 


(40e) 


For the critical problem a^ is an arbitrary constant corresponding to 
the power level and here is normalized to one. At each iteration, a 
search based on successively finer increments to the radius is conducted 
until the criticality condition Eq. (40e) is satisfied. The sequence is 
repeated until the critical radius does not change between successive 
iterations for A n (y). The solution is the same as that of Mitsis with 
exception of the definitions for A (y) , H (v,y) and the inclusion of 
the lQ(R/y) term in Eq. (40b). The reformulation of H(v, y) results in 
a bounded function H (v,y) which leads to a convergent solution. 


III. NUMERICAL METHODS AND RESULTS 
A. Methods 

The various functions and parameters present in the solution were 
calculated by the following methods. The x functions were calculated 
by the method of Shure and Natelson , (id) Eqs. (39a, b,c). The magnitude 
of the discrete eigenvalue was determined by Newton's method. (iiL) Reg- 
ular Bessel functions were calculated with a recursion relationship . (.11) 
The modified Bessel functions were calculated with a series expansion (16) 
followed by an asymptotic series for the larger arguments . 

The calculation was done on the IBM-7094-11 computer. In order to 
achieve single precision accuracy in the results, the calculation was 
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done in double precision. The discrete eigenvalues and the x functions 
agreed with the tabulated values of Kowalska. QZ) The Bessel functions 
agreed with the tabulated values in Abramowitz and Stegun^i^} to at least 
nine significant figures. 

The calculations were done over a quadrature of twenty-four Gauss- 
Legendre points. The variation of the critical radius with the number 
of Gauss points is shown in Table I. This indicates that twenty-four 
Gauss points are sufficient for seven place accuracy. The singular in- 
tegral in Eq. (40c) was evaluated by subtracting \±±) the singularity and 
determining the derivative term with Lagrangian interpolation over all 
points. The critical radius at each iteration for A n (v) was considered 
converged when the value of the integral in Eq. (40e) fell to less than 
10~9. xhe sensitivity of the criticality condition to the variation of 
the radius is shown in figure 2. 


B. Results 

The results presented in Tables II, III, and IV are believed to be 
precise to seven significant figures for the critical radii and to the 
number of figures shown for the neutron density values. Table II lists 
critical radii as a function of c, the mean number of neutron second- 
aries per collision or the infinite medium multiplication factor. Using 
the exact values from the present work as reference values, we see that 
the values of Carlson and Bell^i^ are most accurate for the largest sys- 
tem (c = 1.02) and least accurate for the smallest system (c = 2.0). 
Carlson and Bell used the extrapolated endpoint method for the larger 
systems, R > 1.5. For the smaller systems they interpolated between 
values calculated with the extrapolated endpoint and variational methods. 
The I Tjj method of Hembd^ appears to be particularly effective. The 
reduced accuracy of the I Tq method for large systems is recognized by 
the author. It is interesting to note that for the smallest system 
(c = 2.0), the approximate solutions are more accurate than the 

original reference value. 

The neutron density or total flux was calculated from Eq. (22). 
Table III lists the neutron density relative to the centerline value. 

It is seen that the radial drop-off in the neutron density decreases 
with decrease in the size of the system. Figure 3 shows the neutron 
density distribution for the case c = 2.0 calculated by the method 
presented here and by various order discrete ordinate calculations. The 
discrete ordinate calculations were performed with the TDSN programCl§-) 
using a moment modified quadrature. The neutron density distribution 
from the calculation agrees very closely with the exact distribu- 

tion. However, the excess system multiplication calculated by discrete 
ordinates indicates the increasing importance of the error in the neu- 
tron density distribution as the discrete quadrature order is reduced. 

The first term of the neutron density, Eq. (22), is the asymptotic 
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density arising from the discrete modes. The asymptotic density corre- 
sponds to the diffusion theory solution based on the exact diffusion 
coefficient. Table IV lists the ratio of the asymptotic to total neu- 
tron density as a function of position. The value of this ratio at the 
outer boundary for the case, c.= 1.02- agrees with the corresponding 
value for the Milne Problem'' — ' with c = 1.0. As anticipated, the 
error in the asymptotic density is most severe on the boundary of the 
smallest system. 

In summary, the exact solution for the monoenergetic, bare, in- 
finite cylinder has been presented and compared with approximate solu- 
tions of this problem. It was found results (i) previously used as 
standards were accurate to only three significant figures for small 
systems (for c >1.6). 
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TABLE I. - VARIATION OF CRITICAL RADIUS WITH THE 


NUMBER OF GAUSS-LEGENDRE QUADRATURE POINTS 


Gauss 

c = 2 

c = 1.02 

8 

0.668620870 

9.043254256 

16 

.668612711 

9.043254733 

24 

.668612815 

9.043254733 

40 

.668612853 

9.043254733 


TABLE II. - CRITICAL RADII IN MEAN-FREE PATHS 


9 3 

c Present work Carlson Hendry Hembd 

exact Bell-*- ^3^8 I T 4 

1.02 9.043255 9.0433 — — 9.04458 

1.05 5.411288 5.4118 5.414 5.41152 

1.1 3.577391 3.5783 3.57744 

1.2 2.287209 2.2884 2.28724 

1.4' 1.396979 1.3973 1.39699 

1.6 1.020839 1.0209 1.02085 

1.8 0.807427 0.8067 0.80743 

2.0 0.668613 0.6673 0.670 0.66862 



TABLE III. - NEUTRON DENSITY AS A FUNCTION OF POSITION 


r/R c 

c 



1.02 

1.05 

1.1 

1.2 

1.4 

1.6 

1.8 

2.0 

0.0 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

.25 

.9236 

.9299 

.9361 

.9433 

.9508 

.9550 

.9578 

.9598 

.50 

.7118 

.7340 

.7561 

.7819 

.8093 

.8248 

.8352 

.8426 

.75 

.4124 

.4522 

.4922 

.5399 

.5918 

.6218 

.6421 

.6570 

.85 

.2821 

.3267 

.3718 

.4263 

.4868 

.5223 

.5466 

.5644 

.91 

.2033 

.2492 

.2960 

.3535 

.4181 

.4566 

.4830 

.5026 

.95 

.1502 

.1958 

.2430 

.3016 

.3686 

.4088 

.4366 

.4572 

.98 

.1086 

.1531 

.1999 

.2589 

.3273 

.3688 

.3976 

.4190 

1.0 

.0749 

.1179 

.1641 

.2233 

.2926 

.3351 

.3646 

.3867 


r/R 

c 

TABLE 

IV. - RELATIVE ASYMPTOTIC 

c 

1 TO TOTAL NEUTRON DENSITY 


1.02 

1.05 

1.1 

1.2 

1.4 

1.6 

1.8 

2.0 

0.0 

1.0000 

1.0000 

1.0004 

1.0024 

1.0087 

1.0153 

1.0214 

1.0266 

.25 

1.0000 

1.0001 

1.0006 

1.0030 

1.0102 

1.0174 

1.0239 

1.0294 

.50 

1.0000 

1.0004 

1.0018 

1.0062 

1.0164 

1.0258 

1.0337 

1.0403 

.75 

1.0006 

1.0034 

1.0092 

1.0205 

1.0391 

1.0532 

1.0641 

1.0728 

.85 

1.0033 

1.0110 

1.0226 

1.0406 

1.0652 

1.0821 

1.0946 

1.1043 

.91 

1.0109 

1.0263 

1.0443 

1.0680 

1.0965 

1.1146 

1.1277 

1.1376 

.95 

1.0294 

1.0539 

1.0773 

1.1041 

1.1335 

1.1513 

1.1639 

1.1734 

.98 

1.0775 

1.1084 

1. 1326 

1.1571 

1.1824 

1.1976 

1.2084 

1.2165 

1.0 

1.2318 

1.2337 

1.2368 

1.2424 

1.2522 

1.2601 

1.2664 

1.2716 



R -- 9.04325 R “ 9. 04326 

Figure 2. - Sensitivity of the criticality condition to the 
seventh significant figure in the radius. 
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